
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> rm(list=ls())
> gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 217520 11.7     592000 31.7   592000 31.7
Vcells 319082  2.5    1534954 11.8  1212462  9.3
> require(foreign)
Loading required package: foreign
> require(ineq)
Loading required package: ineq
Warning message:
package 'ineq' was built under R version 3.3.2 
> require(dplyr)
Loading required package: dplyr

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

> library(Hmisc)
Loading required package: lattice
Loading required package: survival
Loading required package: Formula
Loading required package: ggplot2

Attaching package: 'Hmisc'

The following objects are masked from 'package:dplyr':

    combine, src, summarize

The following objects are masked from 'package:base':

    format.pval, round.POSIXt, trunc.POSIXt, units

> require(corrplot)
Loading required package: corrplot
> require(systemfit)
Loading required package: systemfit
Loading required package: Matrix
Loading required package: car
Loading required package: lmtest
Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

> require(cem)
Loading required package: cem
Loading required package: tcltk

How to use CEM? Type vignette("cem")

> require(BayesTree)
Loading required package: BayesTree
> require(stargazer)
Loading required package: stargazer

Please cite as: 

 Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
 R package version 5.2. http://CRAN.R-project.org/package=stargazer 

> require(texreg)
Loading required package: texreg
Version:  1.36.23
Date:     2017-03-03
Author:   Philip Leifeld (University of Glasgow)

Please cite the JSS article in your publications -- see citation("texreg").
Warning message:
package 'texreg' was built under R version 3.3.3 
> require(lavaan)
Loading required package: lavaan
This is lavaan 0.5-22
lavaan is BETA software! Please report any bugs.
Warning message:
package 'lavaan' was built under R version 3.3.2 
> 
> setwd(paste(substr(getwd(),1,regexpr("Dropbox",getwd())[1]-1),"Dropbox/coauthors/AMT_Bisbee_Larson/bisbee_larson_replication_final",sep=""))
> getwd()
[1] "D:/Dropbox/coauthors/AMT_Bisbee_Larson/bisbee_larson_replication_final"
> source("./tools/leg2.R")
> ##############################################################################################
> ##
> ##  File Name:    rels_results.R
> ##
> ##  Input Files:  rels_donate_cem.dta
> ##                final_data_long.dta
> ##  Output Files: rels_donate_cem.pdf             (Figure 8)
> ##                SI_sem_donate.pdf               (SI Figure 4)
> ##
> ##  Purpose: This file summarizes the relationship results for donation behavior in 
> ##            rels_donate_cem.pdf (Figure 8 in the text) as well as applying structural
> ##            equation modeling (SEM) to the data to generate SI_sem_donate.pdf (Figure 4
> ##            in the supporting information).
> ##
> ##############################################################################################
> 
> betas <- read.dta("./relationships/rels_donate_cem.dta")
> 
> ###############################################################################################
> # rels_donate_cem.pdf: Figure 8
> ###############################################################################################
> betas <- betas[complete.cases(betas),]
> betas$diffs <- betas$betaon - betas$betaoff
> betas$ses <- sqrt(betas$seon^2+betas$seoff^2)
> betas$tstats <- betas$diffs/betas$ses
> coefs <- as.data.frame(cbind(as.numeric(betas$betaon),as.numeric(betas$betaoff)))
> coefs$diff <- abs(as.numeric(coefs[,1])-as.numeric(coefs[,2]))
> ses <- as.data.frame(cbind(betas$seon,betas$seoff))
> y.axis <- length(betas$betaon):1
> 
> order1 <- order(betas$diffs)
> coefs1 <- coefs[order1,]
> ses1 <- ses[order1,]
> labels1 <- betas[order1,"names"]
> betas$tie <- sapply(1:length(betas[,"names"]),function(i) strsplit(betas[,"names"]," ")[[i]][length(strsplit(betas[,"names"]," ")[[i]])])
> betas$tie[which(betas$tie=="Weakes")] <- "Weakest"
> 
> betas$names2 <- sapply(1:nrow(betas),function(i) substr(betas$names[i],1,stop = gregexpr(" ",betas$names[i])[[1]][length(gregexpr(" ",betas$names[i])[[1]])]-1))
> 
> betas1 <- betas %>% select(names2,tie,diffs,ses)
> betas1 <- reshape(betas1,timevar="tie",idvar=c("names2"),direction="wide")
> xlim1 <- c(min(betas$diffs-2*betas$ses),max(betas$diffs+2*betas$ses))
> 
> onmin <- coefs[,1]-2*ses[,1]
> offmin <- coefs[,2]-2*ses[,2]
> minmin <- NA
> for(i in 1:length(onmin)) {
+   minmin[i] <- min(onmin[i],offmin[i])
+ }
> 
> order2 <- order(minmin)
> coefs2 <- coefs[order2,]
> ses2 <- ses[order2,]
> labels2 <- betas[order2,"names"]
> ord <- nrow(betas1):1
> 
> # Sprucing up the figure
> betas2 <- betas %>% select(betaon,seon,betaoff,seoff,names2,tie,tstats)
> betas2 <- reshape(betas2,timevar="tie",idvar=c("names2"),direction="wide")
> #betas2[which(betas2$names2 == "Preferred Interaction" | betas2$names2 == "Tie Strength"),which(grepl("beta",colnames(betas2)))] <- betas2 %>% filter(names2 == "Preferred Interaction" | names2 == "Tie Strength") %>% select(which(grepl("beta",colnames(betas2))))*-1
> betas2[,which(grepl("tstat",colnames(betas2)))] <- abs(betas2[,which(grepl("tstat",colnames(betas2)))])
> y.axis <- 1:length(betas2[,1])
> betas2$names2 <- c("Job Search","Contrib. to Entrep.","Garner Conts.","Personal Gain (% $100)","Group Gain (% $100)","Education Hom.","Religious Hom.","Political Hom.","Class Hom.","Tie Strength","Pers. Crisis","Pers. Success","Prof. Crisis","Prof. Success","Pref. Interaction")
> order <- match(c("Political Hom.","Religious Hom.","Class Hom.","Education Hom.",
+            "Contrib. to Entrep.","Garner Conts.","Personal Gain (% $100)","Group Gain (% $100)",
+            "Pers. Success","Pers. Crisis","Pref. Interaction","Tie Strength",
+            "Job Search","Prof. Success","Prof. Crisis"),betas2$names2)
> 
> betas2 <- betas2[rev(order),]
> pdf("./1_Figures/rels_donate_cem.pdf",height=5,width=6)
> nf <- layout(matrix(c(1,2,3,4,5,6,7,7,7,7,7,7),2,6,byrow=T),widths=c(2,1.4,1.4,1.4,1.4,1.4),heights=c(7,.75))
> #layout.show(nf)
> par(mar=c(2,5,3,0))
> plot(rep(0,length(y.axis)), y.axis, type = "n", axes = F, xlab = "", ylab = "",ylim = c(.5,length(y.axis)+.5), main = "")
> axis(2, at = y.axis, label = betas2[,1], las = 1, tick = FALSE, cex.axis =1.3,pos = .9)
> 
> par(mar=c(2,.5,3,.5))
> for(i in c("Strongest","Strong","Weak","Weakest","Clique")) {
+   min <- min(betas2[,paste("betaon.",i,sep="")]-3*betas2[,paste("seon.",i,sep="")],
+              betas2[,paste("betaoff.",i,sep="")]-3*betas2[,paste("seoff.",i,sep="")])
+   max <- max(betas2[,paste("betaon.",i,sep="")]+7*betas2[,paste("seon.",i,sep="")],
+              betas2[,paste("betaoff.",i,sep="")]+7*betas2[,paste("seoff.",i,sep="")])
+   plot(betas2[,paste("betaon.",i,sep="")],y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.2, 
+        xlim=c(min,max),ylim = c(.5,length(y.axis)+.5), main = i)
+   axis(1, at = seq(min,max,(max-min)/10), 
+        labels = sapply(0:10, function(x) round(min+x*((max-min)/10),1)),tick = T,cex.axis = .75, mgp = c(2,.7,0))
+   segments(betas2[,paste("betaon.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis, 
+            betas2[,paste("betaon.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis, lwd =  1,col=rgb(0,0,0,alpha=150,maxColorValue=255))
+   segments(betas2[,paste("betaon.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis-.05, 
+            betas2[,paste("betaon.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis+.05, lwd =  1,col=rgb(0,0,0,alpha=150,maxColorValue=255))
+   segments(betas2[,paste("betaon.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis-.05, 
+            betas2[,paste("betaon.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis+.05, lwd =  1,col=rgb(0,0,0,alpha=150,maxColorValue=255))
+   
+   segments(betas2[,paste("betaoff.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis, 
+            betas2[,paste("betaoff.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis, lwd =  1,col=rgb(0,0,0,alpha=50,maxColorValue=255))
+   segments(betas2[,paste("betaoff.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis-.05, 
+            betas2[,paste("betaoff.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis+.05, lwd =  1,col=rgb(0,0,0,alpha=50,maxColorValue=255))
+   segments(betas2[,paste("betaoff.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis-.05, 
+            betas2[,paste("betaoff.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis+.05, lwd =  1,col=rgb(0,0,0,alpha=50,maxColorValue=255))
+   points(betas2[,paste("betaon.",i,sep="")],y.axis,pch=21,bg=rgb(0,0,0,alpha=150,maxColorValue=255),col=rgb(0,0,0,alpha=150,maxColorValue=255))
+   points(betas2[,paste("betaoff.",i,sep="")],y.axis,pch=21,bg=rgb(0,0,0,alpha=50,maxColorValue=255),col=rgb(0,0,0,alpha=50,maxColorValue=255))
+   stars <- ifelse(round(2*pt(-abs(betas2[,paste("tstats.",i,sep="")]),df = 450),2) < .05,
+                   paste(round(2*pt(-abs(betas2[,paste("tstats.",i,sep="")]),df = 450),2),"*",sep=""),
+                   round(2*pt(-abs(betas2[,paste("tstats.",i,sep="")]),df = 450),2))
+   text(max,y.axis,stars,cex=.8,adj = .9)
+   abline(v=0, lty = 2) # draw dotted line through 0 for reference line for null significance hypothesis testing
+ }
> 
> #add legend (manually) to identify which dots denote model 1 and which denote model 2
> par(mar=c(0,0,0,0))
> plot(rep(0,length(y.axis)), y.axis, type = "n", axes = F, xlab = "", ylab = "",ylim = c(.5,length(y.axis)+.5), main = "")
> legend("center", c("Offline","Online","*  Sig @ 95%","**  Sig @ MC"), 
+        lwd=c(1,1,NA,NA),horiz=T,pch=c(21,21,NA,NA),col=c("grey60","black",NA,NA),
+        pt.bg = c(rgb(0,0,0,alpha=50,maxColorValue=255),rgb(0,0,0,alpha=250,maxColorValue=255),"black"),bg="white",bty='o', cex=.9)
> dev.off()
null device 
          1 
> 
> 
> 
> 
> #################################################################################################
> ##  SI_sem_donate.pdf: Supporting Information Figure 4
> #################################################################################################
> dat <- read.dta("./tools/final_data_clean.dta")
> 
> dims <- c("jobsearch","contribute","garner","personalgain","grpgain","homo_pol","homo_rel","homo_educ",
+           "homo_ls","pers_cris_synth","pers_succ_synth","prof_cris_synth","prof_succ_synth",
+           "seat_synth2_lst","seat_synth2_mst","interaction")
> 
> coefs <- ses <- matrix(NA,nrow = length(dims),ncol = 5)
> i = 1
> for(tie in c("clique","weakest","weak","strong","strongest")) {
+   model <- paste("std_donate_",tie," ~ ",paste(paste("std_",dims,"_",tie,sep=""),collapse=" + "),sep = "")
+   fit <- sem(model,data = dat)
+   coefs[,i] <- fit@ParTable$est[which(grepl("^~$",fit@ParTable$op))]
+   ses[,i] <- fit@ParTable$se[which(grepl("^~$",fit@ParTable$op))]
+   i = i+1
+ }
Found more than one class "Model" in cache; using the first, from namespace 'MatrixModels'
> 
> rownames(coefs) <- rownames(ses) <- dims
> colnames(coefs) <- colnames(ses) <- c("clique","weakest","weak","strong","strongest")
> 
> dat.long <- read.dta("./tools/final_data_long.dta")
> dims <- c("jobsearch","contribute","garner","personalgain","grpgain","homo_pol","homo_rel","homo_educ",
+           "homo_ls","pers_cris_synth","pers_succ_synth","prof_cris_synth","prof_succ_synth",
+           "seat_synth2_lst","seat_synth2_mst","interaction")
> model.long <- paste("std_donate_ ~ ",paste(paste("std_",dims,"_",sep=""),collapse=" + "),sep = "")
> fit.long <- sem(model.long,data = dat.long)
> coefs.long <- fit.long@ParTable$est[which(grepl("^~$",fit.long@ParTable$op))]
> ses.long <- fit.long@ParTable$se[which(grepl("^~$",fit.long@ParTable$op))]
> 
> ords <- order(coefs.long)
> labs <- c("Job Search","Contrib. to Entrep.","Garner Conts.","Personal Gain (% $100)","Group Gain (% $100)","Political Hom.",
+           "Religious Hom.","Education Hom.","Class Hom.","Pers. Crisis","Pers. Success",
+           "Prof. Crisis","Prof. Success","Least in Common","Most in Common","Pref. Interaction")
> 
> pdf("./1_Figures/SI_sem_donate.pdf")
> nf <- layout(matrix(c(1,2,3,4,5,6,7),1,7,byrow=T),widths=c(2,3,1))
> #layout.show(nf)
> par(mar=c(2,5,3,0))
> plot(rep(0,length(coefs.long)), 1:length(coefs.long), 
+      type = "n", axes = F, xlab = "", ylab = "",main = "")
> axis(2, at = 1:length(coefs.long), label = labs[ords], las = 1, tick = FALSE, cex.axis =.9,pos = .9)
> 
> par(mar=c(2,0,3,1))
> pchs <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,19,21)
> cexs <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,.6,1.2)
> cols <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,rgb(0,0,0,.1),rgb(0,0,0,.8))
> ltys <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,1,1)
> bgs <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,NA,"white")
> plot(coefs.long,1:length(coefs.long),type='n',xlab="Coefficient",ylab="",yaxt='n',xlim=c(-.1,.25),main="Full Data")
> segments(coefs.long[ords] - 2*ses.long[ords],1:length(coefs.long),coefs.long[ords] + 2*ses.long[ords],col=cols,lty=ltys)
> points(coefs.long[ords],1:length(coefs.long),pch=pchs,col=cols,cex=cexs,bg=bgs)
> abline(v=0,lty=2)
> par(mar=c(2,0,3,.05))
> ties <- c("Clique","Weakest","Weak","Strong","Strongest")
> # plot(rep(0,length(coefs.long)), 1:length(coefs.long),xlim=c(-.2,.4),
> #      type = "n", xlab = "", yaxt='n',ylab = "",main = "By Tie")
> yadj <- 0
> for(x in 1:5) {
+   pchs <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,19,21)
+   cexs <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,.6,1.2)
+   cols <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,rgb(0,0,0,.1),rgb(0,0,0,.8))
+   ltys <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,1,1)
+   bgs <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,NA,"white")
+   plot(rep(0,length(coefs.long)), 1:length(coefs.long),xlim=c(-.25,.35),
+        type = "n", xlab = "", yaxt='n',ylab = "",main = ties[x],xaxt='n',bty='n')
+   segments(coefs[ords,x] - 2*ses[ords,x],(1:length(coefs[,1]))+yadj,coefs[ords,x] + 2*ses[ords,x],col=cols,lty=ltys)
+   points(coefs[ords,x],(1:length(coefs[,1]))+yadj,pch=pchs,cex=cexs,col = cols,bg=bgs)
+   #yadj = yadj - .15
+   abline(v=0,lty=2)
+   abline(v = .37,col=ifelse(x == 5,"black",rgb(0,0,0,.3)))
+   abline(v = -.27,col=ifelse(x == 1,"black",rgb(0,0,0,.3)))
+   abline(h = 16.6)
+   abline(h = .4)
+ }
> dev.off()
null device 
          1 
> 
> proc.time()
   user  system elapsed 
   3.15    0.57    3.75 
